Supplemental Materials

1 Additional Predictions (Derived Quantities)

1.1 Cardiac Output

Cardiac output is the product of heart rate (hr) and stroke volume (sv).

We will use a parametric bootstrap to make predictions of cardiac output as a function of log10(body mass) and habitat based on our random-effect models for hr and sv. This time, we focus on the difference in predicted cardiac output for (hypothetical) pairs of species of a range of masses that differ only in their habitat.

n_boot <- 1000

We will make sets of many (1000) predictions for each of the species that occur in both the hr and sv data sets.

# note NAMES HAVE TO BE SAME as passed to glmmTMB() when models were fitted (geesh)
mammal_hr <- readRDS('fitted-models/hr-data.RDS')
hr_mod <- readRDS('fitted-models/hr-re-model.RDS')
mammal_sv <- readRDS('fitted-models/sv-data.RDS')
sv_mod <- readRDS('fitted-models/sv-re-model.RDS')

Make dataset for which to make predictions, pairs of hypothetical species that range in size from about 0.1 to 5000 kg.

cardiac_hr <- cardiac_sv <- data.frame(numb = c(1:100),
                                       mass.kg = c(seq(from = 0.1, to = 5000, by = 100),
                                                   seq(from = 0.1, to = 5000, by = 100)))

cardiac_hr <- cardiac_hr %>%
  mutate(log10.body.mass = log10(mass.kg),
         habitat = c(rep('aquatic', 50),
                     rep('land', 50)),
         animal = sample(mammal_hr$animal, size = 100, replace = TRUE),
         order = sample(mammal_hr$order, size = 100, replace = TRUE),
         genus = sample(mammal_hr$genus, size = 100, replace = TRUE),
         species = sample(mammal_hr$species, size = 100, replace = TRUE),
         log10fr = 0) %>%
  dplyr::select(numb, animal, order, genus, species, log10.body.mass, log10fr, habitat)

cardiac_sv <- cardiac_sv %>%
  mutate(log10.body.mass = log10(mass.kg),
         habitat = c(rep('aquatic', 50),
                     rep('land', 50)),
         order = sample(mammal_sv$order, size = 100, replace = TRUE),
         genus = sample(mammal_sv$genus, size = 100, replace = TRUE),
         species = sample(mammal_sv$species, size = 100, replace = TRUE),
         log10sv = 0) %>%
  dplyr::select(numb, log10sv, log10.body.mass, habitat, order, genus, species) 

For each row of the dataset, make a set of 1000 model predictions of cardiac output. These predictions will include uncertainty in the model parameters (intercept and slopes) as well as variation attributed to order/genus/species. (Residual variation is not included). The mass estimates are different in the different datasets; here we use the mass estimate from the dataset to which the model (hr or sv) was originally fitted.

First we define custom functions to make the predictions

predict_hr <- function(.){
  predict(.,
          newdata = cardiac_hr,
          re.form = NA,
          type = 'response')
}

predict_sv <- function(.){
  predict(.,
         newdata = cardiac_sv,
           re.form = NA,
           type = 'response',
         allow.new.levels = TRUE)
  }

Do the bootstrap(s)

if (!file.exists('fitted-models/hr_habitat_boot.RDS')){
hr_boot <- bootMer(hr_mod, 
                   FUN = predict_hr,
                   nsim = n_boot,
                   re.form = NA,
                   type = 'parametric'
                   )

saveRDS(hr_boot, 'fitted-models/hr_habitat_boot.RDS')
}else{
  hr_boot <- readRDS('fitted-models/hr_habitat_boot.RDS')
}

if (!file.exists('fitted-models/sv_habitat_boot.RDS')){
sv_boot <- bootMer(sv_mod, 
                   FUN = predict_sv,
                   nsim = n_boot,
                   re.form = NA,
                   type = 'parametric')
saveRDS(sv_boot, 'fitted-models/sv_habitat_boot.RDS')
}else{
  sv_boot <- readRDS('fitted-models/sv_habitat_boot.RDS')
}

Add boot results to datasets, combine to one dataset, and compute predicted cardiac output (and its median and 95% bootstrap CI):

cardiac_hr <- cardiac_hr %>%
  mutate(boot_pred_hr = unlist(apply(t(hr_boot$t), 1, list), recursive = FALSE),
         across(where(is.factor), as.character))

cardiac_sv <- cardiac_sv %>%
  mutate(boot_pred_sv = unlist(apply(t(sv_boot$t), 1, list), recursive = FALSE),
         across(where(is.factor), as.character)) 
  
cardiac_preds <- inner_join(cardiac_hr, cardiac_sv,
                         by = c("numb", "habitat")) %>%
  mutate(boot_pred_cardiac = map2(boot_pred_hr, boot_pred_sv, function(.x, .y) mapply(prod, 10^.x, 10^.y))) %>%
  mutate(boot_median_cardiac = map_dbl(boot_pred_cardiac, median, na.rm = TRUE),
         boot_lo = map_dbl(boot_pred_cardiac, quantile, probs = 0.025, na.rm = TRUE),
         boot_hi = map_dbl(boot_pred_cardiac, quantile, probs = 0.975, na.rm = TRUE))
# units: ml/stroke * beats/min

plot results

cardiac_mass <- gf_point(boot_median_cardiac ~ log10.body.mass.x, data = cardiac_preds,
         color = ~habitat) %>%
  gf_errorbar(boot_lo + boot_hi ~ log10.body.mass.x, data = cardiac_preds,
              color = ~habitat) %>%
  gf_labs(title = 'Expected Cardiac Output\nfor the species modeled',
          y = 'Predicted Cardiac Output (mL/min)',
          x = 'Body Mass (kg)') %>%
  gf_theme(scale_color_manual(values = my_colors)) 

ggplotly(cardiac_mass, tooltip = '') 

cardiac_mass

cardiac_mass_log <- gf_point(boot_median_cardiac ~ 10^log10.body.mass.x, data = cardiac_preds,
         color = ~habitat,
         text = ~animal) %>%
  gf_errorbar(boot_lo + boot_hi ~ 10^log10.body.mass.x, data = cardiac_preds,
              color = ~habitat) %>%
  gf_lims(x = c(0.001, 10000)) %>%
  gf_refine(scale_x_continuous(trans='log10',
                               labels = scales::label_comma(accuracy = 0.001,
                                                            drop0trailing = TRUE))) %>%
  gf_refine(scale_y_continuous(trans='log10',
                               labels = scales::label_comma(accuracy = 0.001,
                                                            drop0trailing = TRUE))) %>%
  gf_labs(title = 'Expected Cardiac Output\nfor the species modeled',
          y = 'Predicted Cardiac Output (mL/min)',
          x = 'Body Mass (kg)') %>%
  gf_theme(scale_color_manual(values = my_colors)) 
  

ggplotly(cardiac_mass_log, tooltip = 'text')

cardiac_mass_log

1.1.1 Cardiac Output: Terrestrial vs. Aquatic

We’re also interested in judging whether the predicted cardiac output differs between terrestrial and aquatic species. For this, we may better consider the difference between (hypothetical) species that differ only in habitat.

cardiac_diffs <- cardiac_preds %>%
  select(log10.body.mass.x, habitat, boot_pred_hr, boot_pred_sv, boot_pred_cardiac) %>%
  pivot_wider(names_from = 'habitat',
              values_from = c('boot_pred_cardiac', 'boot_pred_sv', 'boot_pred_hr')) %>%
  mutate(boot_diff_cardiac_tma = purrr::map2(boot_pred_cardiac_land, 
                                                boot_pred_cardiac_aquatic,
                                                ~ .x - .y)) %>%
  mutate(boot_median_cardiac_diff = map_dbl(boot_diff_cardiac_tma, median, na.rm = TRUE),
         boot_cardiac_diff_lo = map_dbl(boot_diff_cardiac_tma, quantile, probs = 0.025, na.rm = TRUE),
         boot_cardiac_diff_hi = map_dbl(boot_diff_cardiac_tma, quantile, probs = 0.975, na.rm = TRUE))

require(ggallin)
gf_line(boot_median_cardiac_diff ~ 10^log10.body.mass.x,
        data = cardiac_diffs) %>%
  gf_ribbon(boot_cardiac_diff_lo + boot_cardiac_diff_hi ~ 10^log10.body.mass.x) %>%
  gf_refine(scale_x_continuous(trans='log10',
                               labels = scales::label_comma(accuracy = 0.001,
                                                            drop0trailing = TRUE))) %>%
  # gf_refine(scale_y_continuous(trans=pseudolog10_trans,
  #                              labels = scales::label_comma(accuracy = 0.001,
  #                                                           drop0trailing = TRUE))) %>%
  gf_labs(title = 'Expected Cardiac Output Difference\nTerrestrial minus Aquatic',
          y = 'Cardiac Output Difference (mL/min)',
          x = 'Body Mass (kg)') 

gf_line(boot_median_cardiac_diff ~ 10^log10.body.mass.x,
        data = cardiac_diffs) %>%
  gf_ribbon(boot_cardiac_diff_lo + boot_cardiac_diff_hi ~ 10^log10.body.mass.x) %>%
  gf_refine(scale_x_continuous(trans='log10',
                               labels = scales::label_comma(accuracy = 0.001,
                                                            drop0trailing = TRUE))) %>%
  gf_refine(scale_y_continuous(trans=pseudolog10_trans,
                               breaks = c(-500000, -50000, -5000, -500, -50, 0, 50, 500, 5000, 50000, 500000),
                               labels = scales::label_comma(accuracy = 0.01,
                                                            drop0trailing = TRUE))) %>%
  gf_labs(title = 'Expected Cardiac Output Difference\nTerrestrial minus Aquatic',
          y = 'Cardiac Output Difference (mL/min)',
          x = 'Body Mass (kg)') 

1.2 Total Ventilation

Total ventilation is the product of breathing rate (fr) and tidal volume (vt).

We will use a parametric bootstrap to make predictions of total ventilation as a function of log10(body mass) and habitat based on our random-effect models for fr and vt.

We will make sets of many (1000) predictions for each of the species that occur in both the fr and vt data sets.

# note NAMES HAVE TO BE SAME as passed to glmmTMB() when models were fitted (geesh)
mammal_fr <- readRDS('fitted-models/fr-data.RDS')
fr_mod <- readRDS('fitted-models/fr-re-model.RDS')
mammal_vt <- readRDS('fitted-models/vt-data.RDS')
vt_mod <- readRDS('fitted-models/vt-re-model.RDS')

vent_species <- intersect(pull(mammal_fr, animal), pull(mammal_vt, animal))

The two datasets have 31 species in common.

Make datasets for which to make predictions, containing these species.

vent_fr <- mammal_fr %>%
  dplyr::filter(animal %in% vent_species) %>%
  dplyr::select(animal, order, genus, species, log10.body.mass, log10fr, habitat) %>%
  mutate(across(where(is.factor), factor))

vent_vt <- mammal_vt %>%
  dplyr::filter(animal %in% vent_species) %>%
  dplyr::select(animal, order, genus, species, log10.body.mass, log10vt, habitat) 

For each row of the dataset, make a set of 1000 model predictions of ventilation. These predictions will include uncertainty in the model parameters (intercept and slopes) as well as variation attributed to order/genus/species. (Residual variation is not included). The mass estimates are different in the different datasets; here we use the mass estimate from the dataset to which the model (fr or vt) was originally fitted.

First we define custom functions to make the predictions

predict_fr <- function(.){
  predict(.,
          newdata = vent_fr,
          re.form = NULL,
          type = 'response')
}

predict_vt <- function(.){
  predict(.,
          newdata = vent_vt,
          re.form = NULL,
          type = 'response')
}

Do the bootstrap(s)

if (!file.exists('fitted-models/fr_boot.RDS')){
fr_boot <- bootMer(fr_mod, 
                   FUN = predict_fr,
                   nsim = n_boot,
                   re.form = NULL,
                   type = 'parametric'
                   )

saveRDS(fr_boot, 'fitted-models/fr_boot.RDS')
}else{
  fr_boot <- readRDS('fitted-models/fr_boot.RDS')
}

if (!file.exists('fitted-models/vt_boot.RDS')){
vt_boot <- bootMer(vt_mod, 
                   FUN = predict_vt,
                   nsim = n_boot,
                   re.form = NULL,
                   type = 'parametric')
saveRDS(vt_boot, 'fitted-models/vt_boot.RDS')
}else{
  vt_boot <- readRDS('fitted-models/vt_boot.RDS')
}

Add boot results to datasets, combine to one dataset, and compute predicted ventilation (and its median and 95% bootstrap CI):

vent_fr <- vent_fr %>%
  mutate(boot_pred_fr = unlist(apply(t(fr_boot$t), 1, list), recursive = FALSE))

vent_vt <- vent_vt %>%
  mutate(boot_pred_vt = unlist(apply(t(vt_boot$t), 1, list), recursive = FALSE))
  
vent_preds <- inner_join(vent_fr, vent_vt,
                         by = c('animal', 'order', 'genus', 'species', 'habitat')) %>%
  mutate(boot_pred_vent = map2(boot_pred_fr, boot_pred_vt, function(.x, .y) mapply(prod, 10^.x, 10^.y))) %>%
  mutate(boot_median_vent = map_dbl(boot_pred_vent, median, na.rm = TRUE),
         boot_lo = map_dbl(boot_pred_vent, quantile, probs = 0.025, na.rm = TRUE),
         boot_hi = map_dbl(boot_pred_vent, quantile, probs = 0.975, na.rm = TRUE),
         mean_mass = ((10^log10.body.mass.x + 10^log10.body.mass.y) / 2) )
# units: breaths/minute (?) * mL / breath ---> mL/min

plot results

vent_mass <- gf_point(boot_median_vent ~ mean_mass, data = vent_preds,
         color = ~habitat,
         text = ~animal) %>%
  gf_errorbar(boot_lo + boot_hi ~ mean_mass, data = vent_preds,
              color = ~habitat) %>%
  gf_labs(title = 'Expected Ventilation\nfor the species modeled',
          y = 'Predicted Ventilation (mL/min)',
          x = 'Body Mass (kg)') %>%
  gf_theme(scale_color_manual(values = my_colors)) 

ggplotly(vent_mass, tooltip = 'text')

vent_mass

vent_mass_log <- gf_point(boot_median_vent ~ mean_mass, data = vent_preds,
         color = ~habitat,
         text = ~animal) %>%
  gf_errorbar(boot_lo + boot_hi ~ mean_mass, data = vent_preds,
              color = ~habitat)%>%
  gf_lims(x = c(0.001, 10000)) %>%
  gf_refine(scale_x_log10(#trans = 'log10',
                              labels = scales::label_comma(accuracy = 0.001,
                                                           drop0trailing = TRUE)),
            scale_y_log10(labels = scales::label_comma(accuracy = 0.001,
                                                            drop0trailing = TRUE))) %>%
  gf_labs(title = 'Expected Ventilation\nfor the species modeled',
          y = 'Predicted Ventilation (mL/min)',
          x = 'Body Mass (kg)') %>%
  gf_theme(scale_color_manual(values = my_colors)) 

ggplotly(vent_mass_log, tooltip = 'text')

vent_mass_log